function z = resid3v(X)
global alpha betta depr sigmaC sigmaL phi govexp omega taumax N T
global kstart kistart popweight mistart transfer deductible restrict
global ALPHA lamda tauK0 kmax Delta1 Delta2

if isreal(X)==0
  z(1) = 1e100;
  return
end

n = size(X,1);
z = zeros(n,1);

Delta1 = X(3*T+1);
Delta2 = X(3*T+2:3*T+1+N-1);
lamda = X(3*T+2+N-1:3*T+1+2*(N-1));
transfer = 0;
deductible = 0;

if restrict == 0
  transfer = X(3*T+2:3*T+1+N-1);
  deductible = 0;
  Delta2 = Delta1;
end
if restrict == 2
  transfer = 0;
  deductible = X(3*T+2:3*T+1+N-1);
  Delta2 = - Delta1;
end

if (findstst(kstart-0.2)*findstst(kmax)<0)
kstst =  fzero('findstst',[kstart-0.2 kmax]);%,1e-7,0,0); %find steady state capital
elseif findstst(kstart-0.2)>0
    kstst=kstart-0.2;
else
    kstst=kmax;
end

rstst = betta^(-1)-1+depr; %r from Euler at the steady state
estst = (rstst/alpha)^(1/(1-alpha))*kstst; % express total effective labor from r=F_k
for jj=2:N
ljpart(jj-1)=popweight(jj)*phi(jj)*lamda(jj-1)^(-sigmaC/sigmaL)*(phi(jj)/phi(1))^(1/sigmaL); %this time l1 is pij*phij*lj
end
lstst=estst/(popweight(1)*phi(1)+sum(ljpart)); %express l1 using the expression for total effective labor and l2=f(lamda,l1)
for jj=2:N
    L2a(jj-1) = lamda(jj-1)^(-sigmaC/sigmaL)*(phi(jj)/phi(1))^(1/sigmaL); %last terms in (11) and f_l
end
wstst = (1-alpha)*kstst^alpha*estst^(-alpha); %w=F_e
cstst = 1/(popweight(1)+sum(popweight(2:N).*lamda))*(kstst^alpha*estst^(1-alpha) - depr*kstst - govexp); %c from resource constraint

if size(X,1)>1
    X=X';
end
klast = [kstart X(1:T-1)]; %k_{t-1}
k = X(1:T);
l = X(T+1:2*T);
gamma = X(2*T+1:3*T); %gamma_t

e=[l' (l'*(lamda.^(-sigmaC/sigmaL).*(phi(2:N)/phi(1))'.^(1/sigmaL)))]*(popweight'.*phi);
e=e';

r = alpha*klast.^(alpha-1).*e.^(1-alpha); %r=F_k
w = (1-alpha)*klast.^(alpha).*e.^(-alpha); %w=F_e                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          
Fke = alpha*(1-alpha)*klast.^(alpha-1).*e.^(-alpha);
Fkk = alpha*(alpha-1)*klast.^(alpha-2).*e.^(1-alpha);
c = 1/(popweight(1)+sum(popweight(2:N).*lamda))*(klast.^alpha.*e.^(1-alpha)+ (1-depr)*klast - k - govexp); %resource constraint

cnext = [c(2:end) cstst];
rnext = [r(2:end) rstst];
lnext = [X(T+2:2*T) lstst];
Fkknext = [Fkk(2:end) alpha*(alpha-1)*k(end).^(alpha-2).*estst.^(1-alpha)];

%check bound on tauk
constrtest = c.^(-sigmaC)./cnext.^(-sigmaC) - betta*(1+ (rnext-depr)*(1-taumax));
ind = max(sign(constrtest),0); %ind = 1 if constraint doesn't bind

%period 0
l0 = l(1);
c0 = c(1);
w0 = w(1);
r0 = r(1);
Fke0 = Fke(1);

%FOC for labor at time 0
mu(1) = (omega*l0^sigmaL*(1+sum(ALPHA.*lamda.^(-sigmaC).*phi(2:N)'/phi(1).*L2a)+(Delta1+sum(Delta2.*phi(2:N)'/phi(1).*L2a))*(1+sigmaL))...
    +c0^(-sigmaC)*(Delta1*kistart(1)+sum(Delta2.*kistart(2:N)))*Fke0*(popweight(1)*phi(1)+sum(ljpart))*(1-tauK0))/(w0*(popweight(1)*phi(1)+sum(ljpart)));

%FOC for consumption at time 0
gammanew(1) =(mu(1)*(popweight(1)+sum(popweight(2:N).*lamda))-c0^(-sigmaC)-sum(ALPHA.*lamda.*(lamda*c0).^(-sigmaC))-(Delta1+sum(Delta2.*lamda))*(1-sigmaC)*c0^(-sigmaC)...
    -sigmaC*c0^(-sigmaC-1)*(Delta1*kistart(1)+sum(Delta2.*kistart(2:N)))*(1+(r0-depr)*(1-tauK0)))/(-sigmaC*c0^(-sigmaC-1));
    %Delta1*(mistart(1)-transfer+deductible)+Delta2*(mistart(2)+transfer+deductible));

%periods t>0
lt = l(2:T);
ct = c(2:T);
wt = w(2:T);
rt = r(2:T);
Fket = Fke(2:T);

%FOC for labor at t>0
mu(2:T) = (omega*lt.^sigmaL*(1+sum(ALPHA.*lamda.^(-sigmaC).*phi(2:N)'/phi(1).*L2a)+(Delta1+sum(Delta2.*phi(2:N)'/phi(1).*L2a))*(1+sigmaL))...
    +gamma(1:T-1).*ct.^(-sigmaC).*Fket*(popweight(1)*phi(1)+sum(ljpart))*(1-taumax))./(wt*(popweight(1)*phi(1)+sum(ljpart)));

%FOC for consumption at t>0
gammanew(2:T) = (mu(2:T)*(popweight(1)+sum(popweight(2:N).*lamda))-(1+(Delta1+sum(Delta2.*lamda))*(1-sigmaC))*ct.^(-sigmaC)-((ct'.^(-sigmaC)*(ALPHA.*lamda.^(1-sigmaC)))*ones(N-1,1))'-gamma(1:T-1)*sigmaC.*ct.^(-sigmaC-1).*(1+(rt-depr)*(1-taumax)))./(-sigmaC*ct.^(-sigmaC-1));

%FOC for labor at the steady state
mustst = omega*lstst^sigmaL*(1+sum(ALPHA.*lamda.^(-sigmaC).*phi(2:N)'/phi(1).*L2a)+(Delta1+sum(Delta2.*phi(2:N)'/phi(1).*L2a))*(1+sigmaL))/(wstst*(popweight(1)*phi(1)+sum(ljpart)));

munext = [mu(2:end) mustst];

z(1:T) = gammanew - gamma;
z(T+1:2*T) = ind.*gamma +(1-ind).*constrtest; %either gamma is zero and the constraint is slack (ind=1), or the constraint holds as equalitity
z(2*T+1:3*T) = mu + betta*gammanew.*cnext.^(-sigmaC).*Fkknext*(1-taumax) - betta*munext.*(1-depr+rnext); %FOC for capital

%present value budget constraint for group 1 (at steady state from period T onwards)
z(3*T+1) = sum(betta.^[0:T-1].*(c.^(-sigmaC).*c- omega*l.^sigmaL.*l)) + betta^T/(1-betta)*(cstst^(-sigmaC)*cstst- omega*lstst^sigmaL.*lstst)...
    -c(1)^(-sigmaC)*(kistart(1)*(1+(r(1)-depr)*(1-tauK0))-transfer+deductible);

%present value budget constraint for group j=2:N (at steady state from period T onwards)
for jj=2:N
z(3*T+jj) = sum(betta.^[0:T-1].*(c.^(-sigmaC)*lamda(jj-1).*c- phi(jj)/phi(1)*omega*l.^sigmaL.*L2(lamda(jj-1),l,phi(jj))))...
    +betta^T/(1-betta)*(cstst^(-sigmaC)*lamda(jj-1)*cstst - phi(jj)/phi(1)*omega*lstst^sigmaL.*L2(lamda(jj-1),lstst,phi(jj)))...
    -c(1)^(-sigmaC)*(kistart(jj)*(1+(r(1)-depr)*(1-tauK0))+transfer+deductible);
end

%FOC for lamda
for jj=1:N-1
    dL2dlamda(jj,:) = - l*sigmaC/sigmaL*lamda(jj)^(-sigmaC/sigmaL-1)*(phi(jj+1)/phi(1))^(1/sigmaL); %f_lamda(lamda,l)
    dL2dlamdastst(jj) = - lstst*sigmaC/sigmaL*lamda(jj)^(-sigmaC/sigmaL-1)*(phi(jj+1)/phi(1))^(1/sigmaL);
end

for jj=1:N-1
z(3*T+N+jj) = sum(betta.^[0:T-1].*(ALPHA(jj)*((lamda(jj)*c).^(-sigmaC).*c-omega*L2(lamda(jj),l,phi(jj+1)).^sigmaL.*dL2dlamda(jj,:))...
    +Delta2(jj)*(c.^(-sigmaC).*c-phi(jj+1)/phi(1)*omega*l.^sigmaL.*dL2dlamda(jj,:))-mu*popweight(jj+1).*(c-w*phi(jj+1).*dL2dlamda(jj,:))))...
    -sum(betta.^[1:T-1].*gamma(1:T-1).*ct.^(-sigmaC).*Fket*popweight(jj+1)*phi(jj+1).*dL2dlamda(jj,2:T)*(1-taumax))...
    -c0^(-sigmaC)*(Delta1*kistart(1)+sum(Delta2.*kistart(2:N)))*Fke(1)*popweight(jj+1)*phi(jj+1)*dL2dlamda(jj,1)*(1-tauK0)...
    +betta^T/(1-betta)*(ALPHA(jj)*((lamda(jj)*cstst)^(-sigmaC)*cstst-omega*L2(lamda(jj),lstst,phi(jj+1))^sigmaL*dL2dlamdastst(jj))...
    +Delta2(jj)*(cstst^(-sigmaC)*cstst-phi(jj+1)/phi(1)*omega*lstst^sigmaL*dL2dlamdastst(jj))-mustst*popweight(jj+1)*(cstst-wstst*phi(jj+1)*dL2dlamdastst(jj)));
end

